Il y a 3 principaux types de données spatiales
Ce sont des données mesurées à des points fixes dans l’espace
Données récoltées à des stations météo
geoR, gstat, nlme, glmmTMB, R-INLA
Ce sont des données où c’est la localisation des points dans l’espace qui est aléatoire
L’ensemble des mentions d’occurrence pour la Phragmite commune (Phragmites autralis)
Il s’agit d’un jeu de données concernant le nombre d’espèces de plantes richesse (nb d’espèces) à différents parcelles sur l’île Las Palmas dans les îles Canaries. Cette île est caractérisée par de fortes variations en altitude en raison de la présence de deux volcans sur l’île. Ces données proviennent d’une étude par Irl et al. (2015) et peuvent être téléchargées librement à cette adresse:
https://doi.org/10.5061/dryad.7sk2g/1
Les exemples qui suivent sont fortement inspirés de Zuur et al. (2017).
La base de données légèrement modifiée peut être téléchargée directement à partir de la page github de l’atelier.
d <- read.csv("https://raw.githubusercontent.com/frousseu/introRINLA/master/canary_richness.csv")
head(d)
## Plotname easting northing richness altitude TCI
## 1 AC001 221.1035 3190.789 21 0.7784 1.26
## 2 AC002 221.4035 3192.189 19 0.4473 1.14
## 3 AC003 221.8035 3192.689 25 0.2921 1.28
## 4 AC005 222.1035 3192.889 21 0.0610 1.48
## 5 AC006 222.0035 3192.689 23 0.1424 1.78
## 6 AC007 222.4035 3192.989 22 0.1211 1.86
On a donc un identifiant de parcelle, des coordonnées, la richesse en espèce et deux variables. L’altitude est en kilomètre et TCI est un index de complexité topographique.
En premier, utilisons les fonctions du package sp pour convertir le data.frame en objet spatial. Ensuite, nous téléchargeons l’île en question avec la fonction getData du package raster.
library(sp)
library(raster)
library(scales)
coordinates(d) <- ~easting + northing # convert to spatial
proj4string(d) <- "+proj=utm +zone=28 +ellps=WGS84 +datum=WGS84 +units=km +no_defs" # assign projection
ds <- spTransform(d, CRS("+init=epsg:4326")) # transform to latlon
spa <- getData("GADM", country = "ESP", level = 2) # download Spain as a shapefile
laspalmas <- disaggregate(spa)[225, ] # subset Las Palmas
par(mar = c(3, 3, 0, 0)) # plot locations ans Las Palmas
plot(laspalmas, axes = TRUE)
plot(ds, add = TRUE, col = alpha("darkgreen", 0.5), pch = 16) # alpha adds tranparency to points
On peut également visualiser les localisation avec une carte interactive à l’aide du package mapview.
library(mapview)
mapviewOptions(basemaps = c("Esri.WorldImagery", "Esri.WorldShadedRelief"))
mapview(ds, zcol = "richness")
Avant de faire un modèle, on visualise brièvement les données. Le @data permet d’extraire le data.frame (la table d’attributs) de l’objet d qui est maintenant un objet spatial, plus précisément un SpatialPointsDataFrame.
plot(d@data[, c("richness", "altitude", "TCI")])
Utilisons premièrement un simple GLM avec un effet quadratique sur l’altitude.
fit <- glm(richness ~ altitude + I(altitude^2) + TCI, data = d@data, family = poisson)
summary(fit)
##
## Call:
## glm(formula = richness ~ altitude + I(altitude^2) + TCI, family = poisson,
## data = d@data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8273 -1.3622 -0.2198 0.8298 6.8234
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.89824 0.04299 44.159 <2e-16 ***
## altitude 0.58252 0.06330 9.203 <2e-16 ***
## I(altitude^2) -0.58296 0.03322 -17.550 <2e-16 ***
## TCI 0.67431 0.02908 23.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4493.8 on 889 degrees of freedom
## Residual deviance: 2813.5 on 886 degrees of freedom
## AIC: 6655.6
##
## Number of Fisher Scoring iterations: 5
On peut visualiser le modèle avec le package visreg. L’argument scale permet d’obtenir les prédictions ous l’échelle de la réponse (et non sous l’échelle log).
library(visreg)
par(mfrow = c(1, 2), mar = c(4, 4, 0, 1))
visreg(fit, "altitude", scale = "response")
visreg(fit, "TCI", scale = "response")
Maintenant, plusieurs des parcelles sont très près les unes des autres et on peut se demander si les nombres d’espèces observés sont plus similaires pour les parcelles près les unes des autres. En d’autres mots, reste-t-il de la variation inexpliquée par nos variables explicatives qui pourrait être structurée spatialement? Pour déterminer ceci, on peut faire un variogramme à l’aide des fonctions du package geoR.
Pour ce faire, on peut premièrement extraire les résidus (variation non-expliquée) de notre modèle et les cartographier afin de déterminer visuellement s’il y a des patrons suggérant la présence de corrélation spatiale. Si les petits résidus ou les grands résidus ont tendance à être ensemble, cela nous suggère qu’il y a potentiellement de la corrélation spatiale.
par(mar = c(3, 3, 0, -0))
plot(d, cex = rescale(resid(fit), to = c(0.2, 3)), col = gray(0.5, 0.5), pch = 16,
axes = TRUE)
Il semble y avoir un tel patron, mais cela n’est pas si facile à confirmer visuellement. On peut vérifier cela de façon plus formelle à l’aide d’un variogramme.
Intuitivement, un variogramme représente la variance des différences entre toutes les paires d’observations pour différentes classe de distances.
library(geoR)
v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit))
## variog: computing omnidirectional variogram
plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)",
ylab = "Semivariance")
Si les observations près dans l’espace se ressemble davantage, on devrait s’attendre à ce que les valeurs de semivariance soient plus faibles pour les distances courtes. En fait, une courbe plate suggère qu’il n’y a pas de corrélation spatiale, alors qu’une courbe qui augmente (et qui atteint possiblement un plateau) suggère qu’il y a corrélation.
Nous avons pris les valeurs par défaut de la fonction variog et notamment le nombre de classes de distances qui est relativement faible. Augmentons le nombre de classes afin d’avoir un peu plus de précision sur ce qui se passe à faible distance.
v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit), breaks = seq(0,
20, by = 0.5), max.dist = 20)
## variog: computing omnidirectional variogram
plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)",
ylab = "Semivariance")
On peut voir que la variance est beaucoup moins élevée à de faibles distances. On peut contrôler davantage le variogramme pour inspecter la forme de la relation.
v <- variog(coords = cbind(d$easting, d$northing), data = resid(fit), breaks = seq(0,
10, by = 0.5), max.dist = 10)
## variog: computing omnidirectional variogram
plot(v, main = "Empirical Variogram for Species Richness", type = "b", xlab = "Distance (km)",
ylab = "Semivariance")
Il semble donc y avoir une importante corrélation jusqu’à environ 2km, après quoi les points sont plutôt plats, indiquant qu’au delà de cette distance, la corrélation est faible ou inexistante.
On a donc une corrélation spatiale et il faut donc en tenir compte dans notre analyse! Pour cela, il faut utiliser des méthodes spéciales qui sont disponibles dans les packages geoR et gstat. Cependant, ces packages peuvent être utilisés avec des modèles gaussiens, mais pas avec les GLM (à part geoRglm).
INLA facilite grandement l’inclusion d’effets spatiaux notamment dans les GLMs ou GLMMs auxquels nous sommes habitués.
Différentes valeurs de paramètres
Un variogramme empirique est un variogramme généré à partir de données.
## variog: computing omnidirectional variogram
Les variogrammes théoriques sont les modèles qui sont ajustés aux variogrammes empiriques.
library(gstat)
show.vgms(models = c("Exp", "Gau", "Sph", "Cir", "Mat"), as.groups = TRUE)
Différents types de variogrammes:
Exp: exponentielGau: gaussienSph: sphériqueCir: circulaireMat: MatérnLes différents types de variogrammes sont décrits par des paramètres qui peuvent être ajustés.
show.vgms(kappa.range = c(0.5, 1, 2, 5), models = "Mat", nugget = c(0.1, 0.2,
0.3, 0.4), max = 10, as.groups = TRUE)
Un variogramme est souvent décrit par:
v <- show.vgms(range = 2, nugget = 3, sill = 6, models = "Exp", max = 10, as.groups = TRUE,
ylim = c(0, 11), plot = FALSE)
plot(semivariance ~ distance, data = v[-1, ], type = "l", ylim = c(0, 10), col = "blue",
xaxs = "i", yaxs = "i", xlab = "distance", ylab = "semivariance")
abline(6, 0, lty = 3)
abline(9, 0, lty = 3)
abline(3, 0, lty = 3)
abline(v = 2, lty = 3)
INLA veut dire Integrated Nested Laplace Approximation. Pour une description “accessible” de la théorie derrière INLA, je vous suggère ce lien dans un blogue.
Bayes’ rule
\[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\]
Lorsque appliquée aux modèles
\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]
\(data\): données \(\theta\): paramètres du modèle
SPDE est pour Stochastic Partial Differential Equation
Cette approche se base entre autres sur une discrétisation du Gaussian Random Field (GRF) par un Gaussian Markov Random Field (GMRF)
Utilise une fonction de Matérn pour capturer la dépendance entre les localisations.
library(INLA)
m <- inla(richness ~ altitude + I(altitude^2) + TCI, data = d@data, family = "poisson")
summary(m)
##
## Call:
## c("inla(formula = richness ~ altitude + I(altitude^2) + TCI, family = \"poisson\", ", " data = d@data)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.3899 0.7240 0.3430 3.4568
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.8980 0.0430 1.8138 1.8980 1.9826 1.8978 0
## altitude 0.5828 0.0633 0.4589 0.5826 0.7074 0.5823 0
## I(altitude^2) -0.5831 0.0332 -0.6487 -0.5830 -0.5182 -0.5827 0
## TCI 0.6744 0.0291 0.6169 0.6745 0.7311 0.6747 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 4.078(0.00)
## Number of equivalent replicates : 218.25
##
## Marginal log-Likelihood: -3348.80
m$summary.fixed[, c(1:3, 5)]
## mean sd 0.025quant 0.975quant
## (Intercept) 1.8980388 0.04300921 1.8137909 1.9825871
## altitude 0.5827786 0.06333174 0.4588786 0.7074154
## I(altitude^2) -0.5830944 0.03323408 -0.6486835 -0.5182087
## TCI 0.6743690 0.02910015 0.6168778 0.7311290
cbind(summary(fit)$coef[, 1:2], confint(fit))
## Estimate Std. Error 2.5 % 97.5 %
## (Intercept) 1.8982431 0.04298664 1.8142114 1.9827229
## altitude 0.5825224 0.06329847 0.4589618 0.7070942
## I(altitude^2) -0.5829604 0.03321702 -0.6484524 -0.5182393
## TCI 0.6743125 0.02908278 0.6169035 0.7309101
L’estimation d’un modèle avec un effet spatial avec INLA requiert de passer par plusieurs étapes.
Voici un résumé graphique de ces différentes étapes avec les fonction associées.
La première étape consiste à créer une grille (mesh) qui va être utilisée pour approximer le champs Gaussien. En premier, il faut obtenir les localisations.
locs <- coordinates(d)
head(locs)
## easting northing
## 1 221.1035 3190.789
## 2 221.4035 3192.189
## 3 221.8035 3192.689
## 4 222.1035 3192.889
## 5 222.0035 3192.689
## 6 222.4035 3192.989
Par la suite, on utilise la fonction inla.mesh.2d pour créer la grille. Optionnellement, on peut utiliser une grille non-convexe qui respecte davantage le contour de nos observations.
library(INLA)
hull <- inla.nonconvex.hull(locs, convex = -0.05)
mesh <- inla.mesh.2d(loc = locs, offset = c(1, 5), max.edge = c(1, 3), cutoff = 1,
boundary = hull)
par(mar = c(0, 0, 6, 0))
plot(mesh, asp = 1)
points(locs, col = alpha("black", 0.6), pch = 16, cex = 0.7)
offset
Spécifie l’étendue de l’extension de la grille au delà des observations et l’extension d’une zone tampon au-delà de cette grille. Cette dernière permet de réduire les effets de bordure lors des estimations. L’étendue de cette zone tampon doit être au moins aussi grande que l’étendue de la corrélation spatiale.
max.edge
Permet de spécifier les dimensions maximales des triangles de la grille à l’intérieur et à l’extérieur dans la zone tampon. Plus ces triangles sont petits, plus les approximations sont précises, mais plus les calculs sont longs. En général, il n’est pas nécessaire que les triangles dans la zone tampon soient aussi petits qu’à l’intérieur de la grille.
cutoff
Par défaut, chaque observation sera utilisée comme un coin d’un triangle (vertex). Pour éviter la création de trop nombreux petits triangles, on spécifie une valeur à cutoff au delà de laquelle des points voisins seront ignorés lors de la création des triangles.
boundary
Cet argument permet d’utilier l’étendue non-convexe autour des observations.
Il faut éviter d’avoir des triangles avec des aigles trop aigus. Autrement, les estimations sont de moins bonnes qualités.
Ceci permet de relier les localisations à la grille et d’établir une pondération qui permet d’estimer les valeurs du champs spatial pour chaque localisation. Les localisations situées sur les coins (vertex) seront estimées à partir des valeurs de la grille et les localisations à l’intérieur du triangle seront estimées à partir d’une moyenne pondérée calculée en utilisant les trois coins du triangle dans lequel elles se trouvent.
A <- inla.spde.make.A(mesh, locs)
Dans cet exemple, l’observation en rouge est située à l’intérieur d’un tiangle et les valeurs représentent la pondération qui sera utilisée sur les valeurs de chaque coin pour estimer la valeur de ce point.
par(mar = c(0, 0, 5, 0))
plot(mesh, asp = 1, xlim = c(223, 226), ylim = c(3172, 3176))
points(locs, pch = 16, cex = 2)
i <- 889 # ligne contenant la localisation en rouge
points(locs[i, , drop = FALSE], col = "red", pch = 16, cex = 2)
w <- which(A[i, ] > 0) # pondérations associées au point
points(mesh$loc[w, 1:2], label = 1:nrow(mesh$loc), pch = 1, cex = 4, lwd = 2,
col = "red")
text(mesh$loc[w, 1:2], label = round(A[i, w], 2), adj = c(-0.6, 0.5), font = 2,
col = "red")
Ceci permet de définir les éléments du SPDE et les éléments associés aux caractéristiques du champs spatial.
spde <- inla.spde2.pcmatern(mesh, alpha = 2, prior.range = c(2, 0.5), prior.sigma = c(5,
0.01))
L’argument alpha est pour spécifier une des paramètres de la fonction de Matérn. Ce paramètre doit être fixé entre 0 et 2. Nous devons également définir les priors du champs spatial. L’approche par défaut est d’utiliser des quantiles pour définir les priors sur l’étendue (range) et l’écart-type (standard deviation, sd) associés au champs spatial. Cette façon de faire se base sur l’approche des Penalized-complexity priors développée par Simpson et al. (2017). En résumé, cette approche permet de pénaliser l’étendue de la corrélation spatiale vers l’infini (ce qui réduit la complexité du phénomène) et la variance vers 0 (ce qui réduit également la complexité du phénomène). En pratique, le prior pour l’étendue \(r\) de la corrélation est spécifié avec \(\alpha\) et \(r_{0}\):
\[P(r<r_{0})= \alpha\]
Où \(\alpha\) représente la probabilité que \(r\) soit inférieur à \(r_{0}\). Dans notre cas, prior.range=c(2,0.5) veut dire que nous croyons qu’il y a 50% de chances que le l’étendue de la corrélation spatiale soit supérieure à 2km et donc qu’il y a également 50% des chances q’elle soit inférieure à 2km. Le prior pour la variance du champs spatial, \(\sigma\), est spécifié avec \(\alpha\) et \(\sigma_{0}\):
\[P(\sigma>\sigma_{0})= \alpha\]
Où \(\alpha\) représente la probabilité que \(\sigma\) soit inférieur à \(\sigma_{0}\). Dans notre cas, cela indique qu’on croit qu’il y a 1% des chances que la variance spatiale en richesse soit supérieure à 5. Ici, il faut nous rappeler que nous travaillons sous l’échelle log et que \(e^5\approx148\) espèces.
Ceci permet de définir un index qui sera utile pour spécifier ou récupérer les éléments associés au champs spatial.
spatial.index <- inla.spde.make.index(name = "spatial", n.spde = spde$n.spde)
Le stack est une façon compliquée de fournir les données, les variables et les effets à INLA. Ceci n’est pas essentiel pour des modèles simples, mais ça le devient lorsque les modèles sont plus compliqués ou lorsque l’on veut par exemple générer des prédictions à partir de notre modèle.
d$altitude2 <- d$altitude^2
v <- c("TCI", "altitude", "altitude2")
X <- data.frame(Intercept = 1, d@data[, v])
stack <- inla.stack(data = list(richness = d$richness), A = list(A, 1), effects = list(spatial = spatial.index,
as.list(X)), tag = "est")
On créé la formule décrivant le modèle souhaité. Notez que l’intercept est créé à la mitaine afin de.
model <- richness ~ -1 + Intercept + altitude + altitude2 + TCI + f(spatial,
model = spde)
Finalement, on fait tourner le modèle en appelant la fonction inla et en fournissant les différents arguments. La spécification de compute=TRUE fait en sorte que les distributions postérieures seront calculées pour toutes les observations. La spécification de config=TRUE permettra plus loin de faire des simulations et de générer des valeurs aléatoires à partir de notre modèle.
m <- inla(model, data = inla.stack.data(stack), control.predictor = list(A = inla.stack.A(stack)),
family = "poisson")
summary(m)
##
## Call:
## c("inla(formula = model, family = \"poisson\", data = inla.stack.data(stack), ", " control.predictor = list(A = inla.stack.A(stack)))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 6.8725 23.1234 0.5559 30.5519
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 1.9972 0.1018 1.7965 1.9973 2.1967 1.9977 0
## altitude 0.6313 0.1694 0.2989 0.6311 0.9646 0.6307 0
## altitude2 -0.5258 0.0804 -0.6837 -0.5259 -0.3678 -0.5260 0
## TCI 0.4624 0.0511 0.3617 0.4625 0.5623 0.4628 0
##
## Random effects:
## Name Model
## spatial SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for spatial 2.1257 0.3607 1.4873 2.1031 2.902 2.0639
## Stdev for spatial 0.4277 0.0272 0.3772 0.4267 0.484 0.4245
##
## Expected number of effective parameters(std dev): 216.06(7.583)
## Number of equivalent replicates : 4.119
##
## Marginal log-Likelihood: -2894.36
Le modèle généré (m) est un objet complexe contenant énormément d’information où il n’est pas toujours facile de s’y retrouver. Heureusement, comme nous venons de le voir la fonction summary peut être utilisée pour un sommaire rapide. Si on veut extraire plus d’informations sur notre modèle, la première étape est probablement d’inspecter les différents éléments de la liste formant m.
names(m)
## [1] "names.fixed" "summary.fixed"
## [3] "marginals.fixed" "summary.lincomb"
## [5] "marginals.lincomb" "size.lincomb"
## [7] "summary.lincomb.derived" "marginals.lincomb.derived"
## [9] "size.lincomb.derived" "mlik"
## [11] "cpo" "po"
## [13] "waic" "model.random"
## [15] "summary.random" "marginals.random"
## [17] "size.random" "summary.linear.predictor"
## [19] "marginals.linear.predictor" "summary.fitted.values"
## [21] "marginals.fitted.values" "size.linear.predictor"
## [23] "summary.hyperpar" "marginals.hyperpar"
## [25] "internal.summary.hyperpar" "internal.marginals.hyperpar"
## [27] "offset.linear.predictor" "model.spde2.blc"
## [29] "summary.spde2.blc" "marginals.spde2.blc"
## [31] "size.spde2.blc" "model.spde3.blc"
## [33] "summary.spde3.blc" "marginals.spde3.blc"
## [35] "size.spde3.blc" "logfile"
## [37] "misc" "dic"
## [39] "mode" "neffp"
## [41] "joint.hyper" "nhyper"
## [43] "version" "Q"
## [45] "graph" "ok"
## [47] "cpu.used" "all.hyper"
## [49] ".args" "call"
## [51] "model.matrix"
Les premiers éléments qui devraient nous intéresser sont les distributions postérieures des coefficients associés à nos différentes variables. On peut extraire les quantiles de celles-cià partir du sommaire des effets fixes.
m$summary.fixed
## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.9971786 0.10184005 1.7965115 1.9973494 2.1967014
## altitude 0.6313093 0.16943925 0.2989201 0.6310948 0.9645592
## altitude2 -0.5258307 0.08042505 -0.6837271 -0.5258781 -0.3678250
## TCI 0.4624021 0.05107928 0.3617380 0.4625302 0.5622610
## mode kld
## Intercept 1.9976856 7.381960e-07
## altitude 0.6306687 6.299124e-07
## altitude2 -0.5259636 3.673678e-07
## TCI 0.4627912 2.249087e-06
On peut également représenter graphiquement ces distributions à l’aide des distributions marginales complètes. Comme on peut le voir, tous les coefficients sont relativement loin de zéros.
par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))
invisible(lapply(names(m$marginals.fixed), function(i) {
p <- m$marginals.fixed[[i]]
plot(p[, 1], p[, 2], type = "l", xlab = i, ylab = "density")
abline(v = 0, lty = 3)
}))
Le sommaire des paramètres associés au champs spatial s’accèdent par le sommaire des “hyperparamètres” (hyperparameters). Avec INLA, les paramètres fixes représentent les paramètres de la régression et tous les paramètres de variances ou associés au champ spatial sont considérés comme des hyperparamètres.
m$summary.hyperpar
## mean sd 0.025quant 0.5quant 0.975quant
## Range for spatial 2.1257011 0.36066113 1.4872934 2.1031495 2.9020523
## Stdev for spatial 0.4277356 0.02722437 0.3772038 0.4267148 0.4839872
## mode
## Range for spatial 2.0639373
## Stdev for spatial 0.4245023
On peut également illustrer les distributions postérieures de ces paramètres.
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
invisible(lapply(names(m$marginals.hyperpar), function(i) {
p <- m$marginals.hyperpar[[i]]
plot(p[, 1], p[, 2], type = "l", xlab = i, ylab = "density")
}))
Le champs spatial peut être visualiser de différentes façons. La façon la plus simple est d’utiliser les fonction de INLA pour projeter les valeurs du champs sur une grille. Ceci peut être fait de cette façon.
library(fields)
library(viridisLite)
# https://haakonbakka.bitbucket.io/btopic108.html#92_plotting_the_spatial_mean_field
xlim <- range(d$easting)
ylim <- range(d$northing)
# - Can project from the mesh onto a 300x300 plotting grid
proj <- inla.mesh.projector(mesh, xlim = xlim, ylim = ylim, dims = c(300, 300))
# - Do the projection
mfield <- inla.mesh.project(projector = proj, field = m$summary.random[["spatial"]][["mean"]])
sdfield <- inla.mesh.project(projector = proj, field = m$summary.random[["spatial"]][["sd"]])
par(mfrow = c(1, 2), mar = c(3, 3, 2, 3))
image.plot(list(x = proj$x, y = proj$y, z = mfield), col = viridis(100), asp = 1)
axis(1)
axis(2)
image.plot(list(x = proj$x, y = proj$y, z = sdfield), col = viridis(100), asp = 1)
axis(1)
axis(2)
Les prédictions dans INLA sont générées en soumettant des observations pour lesquelles la variable réponse est NA (dans notre cas, richness=NA). Si on utilise l’argument compute=TRUE, INLA se chargera de calculer des distributions postérieures pour l’ensemble des valeurs avec NA. Cela peut nous permettre par exemple de soumettre des valeurs d’altitude pour étudier son effet sur la richesse.
On génère 50 valeurs d’altitude entre les valeurs minimales et maximales.
n <- 50
x <- seq(min(d$altitude), max(d$altitude), length.out = n)
newX <- data.frame(Intercept = 1, TCI = mean(d$TCI), altitude = x, altitude2 = x^2)
Pour l’exercice, on va assumer que l’effet spatial est constant en prenant un endroit quelconque sur l’île. On répète cette localité autant de fois que le nombre de valeurs d’altitude. Idéalement, il serait plus logique de ne pas tenir compte de l’effet spatial dans les prédictions, mais cela complexifie passablement le code nécessaire.
newlocs <- matrix(c(224, 3175), ncol = 2)[rep(1, n), ]
On associe la localisation fictive à la grille.
A.pred <- inla.spde.make.A(mesh = mesh, loc = newlocs)
On construit une stack pour fournir les valeurs à prédires et on la joint à celle générée plus haut qui contient les données en tant que tel. On spécifie tag=pred ce qui va nous permettre d’identifier les lignes contenant les prédictions.
stack.pred <- inla.stack(data = list(richness = NA), A = list(A.pred, 1), effects = list(spatial = spatial.index,
as.list(newX)), tag = "pred")
stack.full <- inla.stack(stack, stack.pred)
Avec ce tag=pred, on construit un index qui nous permettra de récupérer les lignes contenant les prédictions.
index.pred <- inla.stack.index(stack.full, tag = "pred")$data
Finalement, on fait tourner notre modèle en prenant soin de spécifier compute=TRUE pour que INLA calcule les distributions postérieures. On spécifie également link=1 pour s’assurer que INLA retransforme les valeurs sous l’échelle du nombre d’espèces et non sous l’échelle log.
m <- inla(model, data = inla.stack.data(stack.full), control.predictor = list(A = inla.stack.A(stack.full),
compute = TRUE, link = 1), family = "poisson")
On récupère les valeurs prédites à partir du sommaire des valeurs prédites (summary.fitted.values) en utilisant l’index créé plus haut.
p <- m$summary.fitted.values[index.pred, c("0.025quant", "mean", "0.975quant")]
plot(x, p[, "mean"], ylim = range(unlist(p)), type = "l", xlab = "Altitude en km",
ylab = "Richesse")
lines(x, p[, "0.025quant"], lty = 3)
lines(x, p[, "0.975quant"], lty = 3)
Select the menu on the left to expand / collapse table of contents (TOC) entries. Press button below to collapse all TOC except the top level headings.
If you only want to collapse level 3 headings press this button.